/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package PathModel;

/**
 *
 * @author Alwyn
 */
public class SphericalEarth implements EarthModel {

    private final double R=6371.228; //Radius of earth in km
    private final double DEG_TO_RAD=Math.PI/180.0;

    public double getMeanEarthRadius()
    {
        return R;
    }

    public Location getPointLocation(Location start, double dDistance, double dBearing)
    {
        double dLatitude=start.latitude * DEG_TO_RAD;
        double dLongitude=start.longitude * DEG_TO_RAD;

        double lat2 = Math.asin(Math.sin(dLatitude) * Math.cos(dDistance / R)
                + Math.cos(dLatitude) * Math.sin(dDistance / R) * Math.cos(dBearing));
        double lon2 = dLongitude + Math.atan2(Math.sin(dBearing)
                * Math.sin(dDistance / R) * Math.cos(dLatitude), Math.cos(dDistance / R)
                - Math.sin(dLatitude) * Math.sin(lat2));

        return (new Location(lat2/DEG_TO_RAD, lon2/DEG_TO_RAD));
    }

    public double getDistance(Location start,Location finish)
    {
        double est_dist, longitude1, longitude2, latitude1, latitude2, zeta, delta;

        /* Algorithm must traps very small distances to avoid errors */
        est_dist = Math.pow((finish.latitude - start.latitude), 2)
                + Math.pow((finish.longitude - finish.longitude), 2);
        est_dist = Math.pow(est_dist, 0.5);

        if (est_dist < 0.00001) {
            return (0.0);
        }

        /* convert the points from degrees to radians */
        latitude1 = start.latitude * DEG_TO_RAD;
        latitude2 = finish.latitude * DEG_TO_RAD;
        longitude1 = start.longitude * DEG_TO_RAD;
        longitude2 = finish.longitude * DEG_TO_RAD;

        zeta = Math.cos(longitude2 - longitude1) * Math.cos(latitude1) * Math.cos(latitude2) + Math.sin(latitude1) * Math.sin(latitude2);

        delta = Math.PI / 2 - Math.atan(zeta / Math.sqrt(1 - zeta * zeta));

        return (delta * R);
    }

    public double getInitialBearing(Location start, Location finish)
    {
        double diff, term1, term2;

        double lon1 = start.longitude * DEG_TO_RAD;
        double lat1 = start.latitude * DEG_TO_RAD;
        double lon2 = finish.longitude * DEG_TO_RAD;
        double lat2 = finish.latitude * DEG_TO_RAD;

        diff = lon2 - lon1;
        term1 = Math.sin(diff);
        term2 = (Math.cos(lat1) * Math.tan(lat2)) - (Math.sin(lat1) * Math.cos(diff));
        return (Math.atan2(term1, term2));
    }

    public Location getPathCentre(Location start, Location finish)
    {
        double lat1=start.latitude * DEG_TO_RAD;
        double lat2=finish.latitude * DEG_TO_RAD;
        double lon1=start.longitude * DEG_TO_RAD;
        double dLon=(finish.longitude*DEG_TO_RAD)-lon1;

        double Bx = Math.cos(lat2) * Math.cos(dLon);
        double By = Math.cos(lat2) * Math.sin(dLon);
        double lat3= Math.atan2(Math.sin(lat1)+Math.sin(lat2),
                Math.sqrt( (Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx) + By*By) );

        double lon3= lon1 + Math.atan2(By, Math.cos(lat1) + Bx);

        return new Location(lat3/DEG_TO_RAD,lon3/DEG_TO_RAD);
    }
}
